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ABSTRACT 

Finite difference methods are used to solve the Biot equations for wave 
propagation in a porous medium. The computational domain is a two dimensional 
grid of uniform spacing where truncation of the grid on all sides is accomplished by 
applying homogeneous Dirichlet boundary conditions. The difference method is 
second order in space and time, and is seen to accurately predict phase speeds of the 
primary compressional and shear waves. 
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I. INTRODUCTION 


Mine warfare has resurfaced as a major concern for the United States after the 
Gulf War. The Tripoli and the Aegis cruiser Princeton suffered mine strikes, driving 
home the danger of mines and the need for effective intelligence and mine 
countermeasures. In 1992, the Chief of Naval Operations reorganized the Navy 
headquarters establishing a Director, Expeditionary Warfare (N85). At this time, the 
Navy also established the Program Executive Office for Mine Warfare - reporting directly 
to the Assistant Secretary of the Navy for research, development, and acquisition. [Ref. 1] 
The United Nations estimates that during recent civil and international strife, more 
than 100 million mines have been laid in 62 countries. A United States government 
report estimates that the worldwide total of deployed land mines increases by 500,000 to 
1 million each year. After troops withdraw, land mines and unexploded ordnances 
remai n in the ground, killing and maiming scores of people - farmers, travelers, workers, 
soldiers - every day indiscriminately. The Secretary general of the United Nations 
Boutros Boutros-Ghali now urges stronger political and financial support for efforts to 
neutralize mines worldwide. [Ref. 2] 

In search of a better means for buried mine detection, BBN Systems and 
Technology completed a study experimenting with Rayleigh waves for mine detection in 
1992. The results of their study claimed the probability of detection for anti-tank mines is 
from 74% to 80% and 17 % for anti-personnel mines . The maximum detection range 
reported was found to be 15 feet. [Ref. 3] 



Currently, the Naval Postgraduate School is conducting experiments based on 
Rayleigh wave propagation as mentioned in the BBN report. The goal of the research is 
to improve detection range and the probability of detection of buried objects. 

Finite difference methods applied to Biot theory is intended to enhance the 
experimental research mentioned. Very little work has been done in the past regarding 
Rayleigh wave propagation in porous media. Asymptotic solutions for high and low 
frequencies were found by Deresiewicz in 1962. However, the author has been unable to 
find numerical experiments on Rayleigh wave propagation in porous medium valid for 
intermediate frequency ranges. 

The original goal of this thesis was to produce a finite difference code capable of 
modeling wave propagation in porous media excited by surface tractions. Due to 
numerical instabilities in discretization of the free surface boundary conditions, the 
problem modeled was amended to consider body wave propagation in porous media. A 
two dimensional solution to the problem is discretized using an evenly spaced grid for the 
"X" and "Y" coordinates, whereby the grid is truncated by an application of homogeneous 
Dirichlet conditions "relatively" far from the point of application of an incident 
displacement. Because the spatial discretization is uniform, and center differencing is 
used throughout, the numerical scheme is second order in space and time. The code is 


written in C++. 



II. THEORY 


Biot theory is one of two popular theoretical methods to model wave propagation 
in porous media, the other theory assumes the media is viscoelastic. Unlike viscoelastic 
theory, which is based more upon elasticity theory, Biot theory predicts the existence of 
three kinds of body waves, two dilatational and one rotational. The first type of 
dilatational wave involves motions of the skeletal frame and the interstitial fluid which 
are nearly in phase and for which the attenuation owing to viscous losses is relatively 
small. The second type of dilatational wave is such that the frame and fluid components 
are moving largely out of phase. Biot derived two coupled linear differential equations 
governing the relationship between the displacements of the fluid and the frame. [Ref. 4] 

A. WAVE EQUATIONS 

The Biot equations are [Ref. 4]: 

p, 1 0 + P12 = (P - P)V(V • u) + M-V 2 u + QV(V • U) + bF(co)(f - f) (1) 

p 120 + P220 = QV(V • u) + RV(V • U) - bF(co)(f -f) (2) 




u = the displacement vector of the solid material; 

U = the displacement vector of the fluid; 

P = the ratio of the volume of the pores to the total volume of the media; 
|i = the shear modulus of the solid; 

Kb = the bulk modulus of the free-draining frame; 

Kf = the bulk modulus of the fluid; 



K r = the bulk modulus of the solid; 
pr = the mass density of the pore fluid; 
p r = the density of solid material; 
r| = fluid viscosity; 

K = permeability; 
a = added mass term; 

F(co)= frequency dependent function to account for high frequency attenuation; 
such that, 

Q = [|3( 1 - (})K r - pK b M 1 - p + - £] 

R = [p 2 K r ] + [l-p + p£-£] 
b = [riP 2 ] + [k] 

P = [(1 - P)(( 1 - P) - ff)Kr + P^l + [1 - P + Pg - + f H 

P12 = Prp(l - a) 

pll = (1 — P)p r — Pl2 
p22 = PfP Pl2 

Mathematical notation : The bold letters denote vectors and the subscripts denote 
derivatives with respect to the superscripted variable. 



V 2 u = 


Uxx+I 



The u and U terms in equations (1) and (2) can be isolated as follows: 

Uii(p22 - fit) = QV(V • u) + RV(V • U) - bF(<o)(U, - u t ) - {fit 

[(P - n)V( V • u) + |iV 2 u + QV(V • U) + bF(<o)(U, - u,)]}; (3) 

u„(p ,2 - = QV(V • u) + RV(V . U) - bF(w)(U, - u t ) - {ff 

[(P - p) V(V • u) + |iV 2 u + QV(V • U) + bF(eo)(U, - u,)]}; (4) 

Rewriting equations (3) and (4) in scalar form: 

U„(p22 - fit) = Q(u*x + v yx ) + R(U XX + V yx ) - bF(©)(U, - u t ) - {fit 

[(P - p)(u xx + v yx ) + (i(u xx + u yy ) + Q(U XX + V yx ) + bF(ro)(U, — u,)]} (5) 

V„(p22 - fit) = Q(u xy + V yy ) + R(U xy + V yy ) - bF(co)( V, - V,) - {fit 


[(P-W(U X y + V yy ) + M.(V, 


,) + Q(U xy + V yy ) + bF(co)(V, - V,)]} 


( 6 ) 



U«(p 12 - = Q(Uxx + Vyx) + R(Uxx + Vyx) - bF(«))(U, - u.) - 

[(P - |X)(u X x + V yx ) + |i(Uxx + u yy ) + Q(Uxx + Vy X ) + bF(a))(U t - u,)]} (7) 

V,,(p 12 — = Q(Uxy + Vyy) + R(Uxy + Vyy) - bF^XV, - V.) - 

[(P - P)(Uxy + Vyy) + P(V XX + Vyy) + Q(U X y + Vyy) + bF(fi))(V t - V,)]} (8) 

B. BOUNDARY CONDITIONS 

A rectangular grid is used in the X and Y direction. Displacements u and U are 
set to zero at all sides of the numerical domain. This creates reflections as body waves 
strike any of the sides. Originally, on one of the surfaces, a traction free condition was 
applied which is necessary for the propagation of surface waves and for proper modeling 
of a porous half space. On this surface, both the stresses and the fluid pressure vanish, 
resulting in the following boundary conditions: 


CTyy I y=0 = (P - 2N)(Ux) + PVy + Q(U X + Vy)l = f(X, t)j (9) 

Oxyly=0=N(Uy+Vx)ly=0 = 0; (10) 

S I y=0 = Q(Ux + Vy) + R(Ux + Vy)l y=0 = -P X f(X, t); (11) 

Uy + Vx=0; (12) 



where o is the stress tensor and S is related to the fluid pressure (-pP), and f(x,t) is an 
applied normal force on the surface of the half space. Equations (10) and (12) are the 

conditions that on the free surface, shear stress is zero for both the solid and the fluid. 

C. SOURCE FUNCTION 

The source function used for the half space problem was: 

F(x, t) = Ae-k-^V'-^ 2 ; (13) 

The source function is Gaussian in shape along the x grid on the surface. x 0 is the 
midpoint of the x grid, while t 0 is the start-up time for the numerical calculations. 'A' is a 
constant set equal to one pascal. 

When it was determined that the discretization for the free surface conditions were 
unstable, the following two initial states of the displacements were employed to check the 
accuracy and stability of the code at interior points of the domain. 

(14) 

(15) 

(16) 

(17) 

7 




Equation (14) and (15) were used to check the compressional waves. They can be 
seen as an initial expansion of the solid portion of the porous medium with two 
dimensional gaussian shape. Equation (16) and (17) are the cross product of equations 
(14) and (15) with an out of the page vector (k). They represent an infinitesimal rotation 
of the solid elements of the two dimensional gaussian distribution of magnitudes, and are 
used to check the shear wave of the porous medium. Both i 0 and j 0 are set at the center 
of the grid, and gamma is set equal to 0.1. 
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III. PROGRAM ALGORITHM 


The program is written in C++. The spatial grid for both the X and the Y axis are 
chosen to be 0.1 meter, while the time interval. At is 0.1 second. The time interval meets 
the stability requirements to be stated shortly. The Biot equations are scaled prior to the 
running of the program. Three different sets of Biot parameters were used. Table 1 
shows the inputs which remain unchanged for all three parameters sets used. 



Omega (fo) 

2000 rad/s 

F(w) 

1 

Grid spacing (h) 

0.1 

Time interval (delta tl 

0.1 


Table 1. Input Parameters 


A. FINITE DIFFERENCE SCHEME 

In all the finite difference approximations, center differences are used. 
Nodal points are found at ( Xi,yj ) where: 

Xi = iAx; (the horizontal component) 

yj = jAy; (the vertical component) 

and the discrete times at which the displacements are found are t n = nAt; 


9 




The term u"j therefore refers to u(xj,yj,t n ) m , (See the diagram below). 



Some of the difference formula's used are: 


u, = i[< 1 - u"; 1 ] + 0(Ar 2 ); 


«« = ^[«i+w - 2u h + <-1J + 0(/i 2 ); 

“jy = — w i+i j-i + u i-[ j-i ~ u i-ij+i ] + 0{h 2 ), 


The "O" terms are neglected leading to a second order truncation errors throughout the 
code. 


10 



















When equations (5)-(8) are differenced in time and space, we may then isolate the 


unknown displacements at time level 'n+1' as follows: 


+ CUy‘ =E 

(18) 

-Dul"-' = F 

(19) 

+ Cv^'=G 

(20) 

-Dvy'=H 

(21) 


The goal is to solve for each grid point's value at the next time level (n+1), to update (or 
shift time levels), and march forward in steps of At. 

Solving for the 'n+1' terms from equations (18) through (21) we have: 


_ DE+CF 
AD+BC 


( 22 ) 


(23) 


Vff 1 


DG+CH 

AD+BC 


(24) 


BG-AH 

AD+BC 


(25) 




The variables A, B, C, D, E, F, G, H in their respective finite difference form are as 
follows: 


A = bF(cD)^ + £fbF(a))^; 

(26) 

B = ^r(p22 -fit j + bF(<n)^ + £fbF(o>)^; 

(27) 

c = ^ (p 12 -Tfir ) - bF (“)lk - fr|bF((0)^; 

(28) 

D = bF((o)^ + ^bF(Q))^; 

(29) 


E = Q{^( u "+i j ~ 2u "j + j)+ 

4h2 ^ v i+l,j+l -v i-lj+l +v i-l,j-l ~ v i+lJ-l j} + 

R{^(u^j-2Ur j +U". 1J )+ 

^(v" +lJ+ , - v?_ 1J+1 +v:_ 1 , H -vr +1J . l ]>- 

bF( ( o)^{ur,- 1 -U^ 1 }- 

( p - w{^(ur + .,j - 2<j+u^j j+ 
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4h 2[ V i+lJ+I v i-lj+l + v i-lj-l v i+lj-lJ} 

-4u"j + u"_ij + u|j +1 +U? J _ 1 ]}- 

ffQ{^(u n WJ -2U? J + U" W)t )+ 

j^[ v M.,j +1 - v !’-ij + 1 + v "-ij-i -v" + ij_, j } - 

0bF((o)^{u|7' -Ufj 1 }- 

- £ f^ i )( u "7' - 2u|7' j; (30) 

F = Q{^(ul’ +IJ -2ui; j + u| , _ 1 , j )+ 

^r( v i+lj+l — v i-lJ+l +v i-lj-l —v i+lj-l)} + 

R{^(ur +1 j-2U|' J +ur_ 1J )+ 
ir(vr + , J+1 -v^ 1(j+ , + v?_ liH + vr +M _,)}- 


13 



pir( p - M-Hp-( u "+i,j - 2u-j+u?_,j )+ 

4h^ ^ V i+lj+l — v i-lj+l + v i-lj-l —v i+lj-l^} — 

^(“WJ _4u !j + “E-IJ +u" J+1 + uPj., )}- 

^fQ{^(ur +IJ -2ur J+ uiL 1J )+ 

^(vr +1J+1 -v"_ u+I +v"_ 1J _ 1 -vr +1J _,)}- 

< 3i > 

G = Q{^i(ur + , J+1 - u"_,j +1 + uf.,J_. - uf +1J _j )+ 

- 2vrj+v"j_, )>+ 

R ^( u ^u +i -ur-., j+ .+ur_ 1J _ 1 -uiv 1J _ 1 )+ 

^( v r J+ >-2vi’ +v^]}- 
bF ( to )(i)( v S7 , - v !j 1 )- 
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u i-l j+1 u i-lj-l 


P'( V U+ | _ 2vfj + v" j _ 1 ^ }— 

M.{p( v ? +1 j - 4v?j + v"_ uj + v? J+1 + v[j_,) }- 

ftQ{^(ur +1J+1 - ur_ 1 J+ 1 +ur_ 10 _, - ur +IJ _,)+ 

^ +1 -2V" j+ V|V_ 1 )}- 

^(p | 2- £ f^)( v u'-2v|; j ); (32) 

h = ~ u i’-ij+i + u T-i j-i - u !+ij-i )+ 


h 2 ( v i,j+l 2v i,j + v i,j-l ) } + 

R{^(uf tIJ+1 -u|L lJ+1 +ur-.,j-. -u? +1<i _ I )+ 


4 V" i+ , -2V":+V",_ 
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^KikX^'-vu 1 )- 


pif(P - W{^(u? + U + 1 - Uwj +1 + Clj-l - U " + 1J-1 )+ 
tf( V £i+l - 2v"j + vfj_! j}- 

-4vrj + vjL u +vf J+1 + V^_, )}- 

^Q{^(ur + . J+ , -ur_ 1J+l +ur_ ld _, -u^,)+ 
p( V ^> -2Vlj^+V"y-, )}- 
pif b F(“)(ik)(v"j 71 - v !j 1 )- 

^(p22-S)(vl‘J , -2Vj' J ); (33) 


Second order finite difference forms of the boundary conditions equations (9) 
through (12) are given by: 


(P - 2N){^(u? +1J -uJL,j)} + P{i(vr j+1 - v^-i)}+ 
Q^UH.j-Ur-ij+V^, ~ V"j_,)} = f"; (34) 
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} = 0; 


(35) 


Q{a( U WJ“ U MJ + V ?tH _V U-l) + 

R {a( U i+ I .j-U"-ij + V" j+l -V^,) =-pfT; (36) 

^Ug+i _ U"j_i + V" + ij j =0; (37) 

Equations (34) through (37) give the means to calculate the 'j-1' terms exterior to 
the boundary surface (see Figure (2)), which are in turn substituted into equations (30 
through (33) to find the values of u and U along the free surface boundary. 



B. SCALING 

Only two parameters need to be scaled in the Biot equations, time and length. The 
equations for the displacements are linear, therefore their scaling is determined by the 
magnitude of the applied force on the free surface or by the initial displacements 
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assigned. The time scale chosen is omega in rad/sec, while the length scale used is L=1 


meter. These scales lead to the following relationships: 

x=Lx s'- i=W 

ts = ®t; t = fc; | = <- (39) 

In addition to scaling, a constant term of ^ is multiplied across each equation. 
The result is a constant multiplier to all terms involving second derivatives with 
respect to time, ^ to all first derivatives in time, and £ to all second derivatives in x or 

y- 

C. STABILITY 

The system of equations (18) to (21) was found to be stable provided that 


At<h,(v? c+ Vl) 2 ; (40) 

Vpc = . ( 4i) 




(42) 


K c = bulk modulus of the saturated medium 
= pK r + (l-(i)K r ; 
p = Pp f +(l-p)p r ; 

(i= shear modulus of the frame; 

h= the grid interval in both x and the y direction; 

At= the time step; 



V Pl =the saturated velocity of the fast compressional wave; 

V s =the shear wave velocity; 

Equation (40) is the stability condition used in finite-difference algorithms 
involving wave propagation for elastic media. [Ref. 5] For porous medium, equation 
(40) is also valid because the saturated velocity of the fast compressional velocity in the 
porous elastic solid is normally much less than the compressional velocity of the purely 
elastic solid. 

After scaling, the stability equation (40) becomes: 



D. PARAMETERS 

The parameters used are shown in the table below. 
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The Stoll and Kan parameters and the ARL:UT parameters are taken from [Ref. 6]. 
ARL-.UT parameters are believed to be the better parameters to use based on the 
experiments performed by Chotiros, ARL University of Texas at Austin. The modified 
ARL:UT is the same as ARL:UT parameters except the fluid viscosity is set to zero, 
eliminating any damping in the system, and smaller grain density and lower porosity 
ratios are used. The wave speeds calculated using the above parameters are shown in 
Table 3. 


r 




Stoll & Kan ARL:UT 

ARL:UT mod 

Fast P wa\« speed 

1485.14m/s 1694.8m/s 

1676.7 m/s 

Slow P wave speed 

56.21 m/s 192.83 m/s 

1066.8 m/s 

Sheer speed 

119.24m/s 114.7 m/s 

114.07 m/s 

Table 3. Wave Speeds 


k_ 




Equations used to calculate the compressional and shear wave speeds are [Ref. 7]: 
Pfaa = ,/Real((2AH)/(B-jf-Del)) ; 


P s iow = VReal((2AH)/(B - if + Del)) ; 


S = J Real(p(( P 22/p) — jf)/(p(C -jf))) ; 

where, 

H = P + 2Q + R; 

A = (PR-Q 2 )/H 2 ; 

B = (pnR + P22P-2p,2Q)/(pH); 

C = (pnp22-p 2 2 )/p 2 ; 

Del = 7(B-jf) 2 -4A(C-jf) ; 
f = |lp/(Kp(0). 
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E. DISCUSSION 


The source was applied to the finite difference model only once at the initial time 
level in an attempt to model an instantaneous point force surface traction on a porous 
half-space. Because of the porous parameters used, the shear speed was small compared 
to the fast compressional wave. The large difference in speeds is assumed to be the cause 
of a numerical instability at the free surface boundary. An instability of this form was 
reported in the work of Ilan and Loewenthal in 1976 [Ref. 8]. To test this hypothesis the 
value of the shear modulus was increased to within the reported stability region of 
center-differencing of the free surface conditions, but still the code was unstable. 
Alternative differencing of the free surface conditions proved futile, so a decision was 
made to drop the free surface altogether and apply homogeneous Dirichlet Boundary 
Condition on all sides of the numerical domain. This numerical problem identifies an 
area of future research. 

The Stoll & Kan parameters demonstrate short propagation distances due to the 
relatively large frequency and high attenuation of the media. (Figure (3)) The ARL:UT 
which neglects fluid viscosity gives a much better contrast for the fast compressional 
wave propagation as reflections off the sides can also be observed. (Figures (4) and (5)) 
The fast and the slow compressional wave are visible using the ARL:UT mod parameters. 
(Figures (6) and (7)) From Figure 7, the fast and slow compressional wave speeds can be 
calculated using equation (44). 
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co = 2000 rad/sec; 


(44) 


L=1 m. 

For the fast compressional wave speed, it traveled 50 grids in about 60 time level giving 
us a speed of 1666 m/s. The slow compressional wave speed traveled 50 grids in 100 time 
level giving us a speed of 1000 m/s. 

To check the accuracy and stability of the code at interior points of the domain, 
equations (14) through (17) were used and the results shown in Figure (8) and (9). In 
Figure (8), the cylindrical spreading show a propagation of approximately 80 grids in 100 
time level for a 1600 m/s fast compressional wave speed. In Figure (9), no shear 
propagation in the fluid is seen as expected for fluids do not support shear propagation. 
For the solid, shear speed propagated approximately 11 grids in 200 time level giving us a 
shear wave speed of 110 m/s. It is extremely difficult to distinguish the slow wave 
propagation from that of the fast wave in Figure (8). However, due to the fact that the 
frame and fluid components of the slow wave are moving largely out of phase and the fast 
wave moves nearly in phase, a plot of relative speed between the frame and the fluid 
would make the slow wave more visible. (Figure (10)). In Figure (10), the slow wave is 
seen to propagate to about 50 grids in 100 time levels for a slow speed of 1000 m/s. 

All of the figures shown and numbers calculated have propagation speeds that 
match well with the theory shown in Table (3). 
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y grid 



30 






IV. CONCLUSION 


The propagation of the fast compressional waves can be clearly seen in all three 
parameter sets used. The program behaved as expected concerning the body waves and 
their speeds matched well with theory. The different parameters used have shown the 
sensitivity of Biot parameters in wave propagation. 

One area that should be pursued is to determine a stable free surface boundary 
condition difference operator so that Rayleigh wave propagation may be modeled. 
Another improvement on the current code would be to include absorbing boundary 
conditions on the sides and the bottom, to better model a semi-infinite domain. 
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APPENDIX PROGRAM CODE 




/* Thesis Program for Wave Propagation */ 

/* Programmer: Jonah Shen */ 

/* Date: June 1995 */ 

/* Language: C, beta, scaled */ 


/***************************************************/ 

#include <stdio.h> 

#include <math.h> 


/* Initialization of variables and constants */ 


float beta = 0.47; 

/* beta is the ratio of the vol of */ 

/* the pore to total vol of element */ 

float Mu = 2.6 le7; 

/* shear modulaus. Pa */ 

float rhor = 2650.0; 

/* grain density, Kg/m A 3 */ 

float rhof = 1000.0; 

/* fluid density, Kg/m A 3 */ 

float Kb = 4.36e7; 

/* bulk modulus of frame. Pa */ 

float Kr=3.6el0; 

/* bulk modulus of solid, Pa */ 

float Kf = 2.0e9; 

/* bulk modulus of fluid. Pa */ 

float N = le-3; 

/* liquid viscosity, kg/m s */ 

float K = le-10; 

/* liquid permeablility, m A 2 */ 

float delT = 0.1; 

/* delta t */ 

float Fw = 1.0; 

/* forcing function */ 

float h = 0.1; 
float alpha = 3.0; 
float w = 2000.0; 

/* step size */ 


float a,Q,R,P,b,rhol2,rho 11 ,rho22,A A,BB,CC,DD,EE,FF,GG,HH,II,JJ,KK,LL; 
float MM,NN,y,f; 

float Aa,Bb,Cc,Dd,Ee,Ff,Gg,Hh,Ii,Jj,Kk,Ll,Mm,Nn,Oo,Pp,Qq,Rr,ans; 
float aA,bB,cC,dD,eE,fF,gG,hH,iIjJ,kK,lL,mM,nN,oO,pP,qQ,rR; 
int n,i,j,ind,time,tt; 


*/ 

*/ 

*/ 


/************************************************* 

/* Defining Vectors for finite difference 

/* Lower case is displacement of points in the skeletal frame 

/* Upper case is dispalcement of fluid 

/* Each letter will be follow by 3 numbers, the first number 

/* represent time second represent x direction, third represent 
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/* y direction */ 

/**********************************************************************/ 

float u[3][51][51]; 
float v[3][51][51]; 
float U[3][51][51]; 
float V[3][51][51]; 
float sol[300][51]; 
float flu[300][51]; 

/* Declare Function Prototype */ 

float FEE(float Aa.float Bb,float Cc.float Dd,float Ee.float Ff,float Gg, 
float Hh,float Ii,float Jj,float Kk.float Ll.float Mm,float Nn, 
float Oo,float Pp,float Qq.float Rr); 
float FFF(float Aa.float Bb,float Cc.float Dd.float Ee.float Ff.float Gg, 
float Hh.float Ii.float Jj .float Kk.float Ll.float Mm.float Nn, 
float Oo,float Pp.float Qq.float Rr); 

float FGG(float aA.float bB,float cC.float dD.float eE.float fF,float gG, 
float hH,float il,float jj,float kK,float lL.float mM.float nN, 
float oO,float pP.float qQ,float rR); 

float FHH(float aA.float bB,float cC.float dD.float eE.float fF.float gG, 
float hH.float ii.float jJ.float kK.float lL.float mM.float nN, 
float oO,float pP.float qQ.float rR); 


main() 

( 

a = (1 -beta+(be ta*Kr/Kf)-( Kb/Kr)); 

Q = (beta*( 1 -beta)*Kr-beta*Kb)/a; 

R = beta* beta* Kr/a; 

P = ((l-beta)*(( l-beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a + 4*Mu/3; 

b = N*beta*beta/K; 

rhol2 = rhof*beta*(l-alpha); 

rholl = (l-beta)*rhor-rhol2; 

rho22 = rhof*beta-rhol2; 

AA = b*Fw*w/(2*delT*Mu) + rho22*b*Fw*w/(rhol2*2*delT*Mu); 

BB = (w*w/(delT*delT*Mu))*(rho22-(rhol2*rhol2/rholl))+b*Fw*w/(2*delT*Mu)+ 
b*Fw*w*rhol2/(2*delT*Mu*rhol 1); 

CC = (w*w/(delT*delT*Mu))*(rhol2-(rho22*rhol l/rhol2))-b*Fw*w/(2*delT*Mu)- 
b*Fw*w*rho22/(2*delT*Mu*rhol2); 

DD = b*Fw*w/(2*delT*Mu)+rhol2*b*Fw*w/(rholl*2*delT*Mu); 

/* Initialization of the arrays */ 
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for (n=0;n<=2;n++) 

{ for (i=0;i<=50;i++) 

{ for (j=0;j<=50;j++) 

{ u[n][i][j] = 0.0; 

v[n][i][j] = 0.0; 
U[n] [i][j] = 0.O; 
V[n][i][j] = 0.0; 


/* Loops */ 

time=200; 
ind = 0; 


/* Boundary Conditions */ 
for (i=l;i<=49;i++) 

{ j=i; 

y=1.0; 

f=exp(-(i-25)*(i-25)); 

U[nl[il[j-l]=U[n][i][j+l]+V[n][i+l][j]-V[n][i-l]|j]; 

u[nJ[i][j-l]=u[n][i][j+l]-v[n][i-l]|j]+v[n][i+l][j]; 

II=-(2*h*f*y)+((P-2*Mu)*(u[n][i+l]|jl-u[n][i-nDl)+P*v[ n ][i](j+l] 
+Q*( V[n][i] [j+ H-U[n] [i-11 [j]+U[n] [i+1 ] [j])); 
JJ=(beta*2*h*f*y)+Q*(v[ n ][i][j+l]-u[n][i-l]|j]+u[ n ][i+l]Ij])+ 
R*(V[n][i][j+l]-U[n][i-l][j]+U[n][i+l]|j]); 

v[n][i][j-l]=(IFR - Q*JJ)/(P*R - Q*Q); 

V[n][i][j-1J=(Q*IT - P*JJ)/(Q*Q - P*R); 


while (ind < time) 

{ 

for (i=l;i<=49;i++) 

{ for (j=l;j<=49;j++) 

{ 

Aa=u[n][i+l][j]; 
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Bb=u[n][i][j]; 

Cc=u[n][i-l][j]; 

Dd=u [n-1 ] [i] [j]; 

Ee=u[n][i]U+lJ; 

Ff=u[n][i][j-1]; 

Gg=v[n][i-l][j+l]; 

Hh=v[n][i-l][j-l]; 

Ii=v[n][i+l][j+l]; 

Jj=v[n][i+l][j-l]; 

Kk=U[n][i+l][j]; 

Ll=U[n][i][j]; 

Mm=U[n][i-l][j]; 

Nn=U[n-l][i][j]; 

Oo=V[n][i+l][j+l]; 

Pp=V[n][i-l][j+l]; 

Qq=V [n] [i-1 ] [j-1 ]; 

Rr=V[n][i+l][j-l]; 

aA=u[n][i+l][j+l]; 

bB=u[n][i-l][j+l]; 

cC=u[n] [i-1 ] [j-1 ]; 

dD=u[n][i+l][j-l]; 

eE=v[n][i][j+l]; 

fF=v[n][i][j]; 

gG=v[n][i][j-l ]; 

hH=v[n-l][i][j]; 

il=v[n][i+l][j]; 

jJ=v[n][i-l][j]; 

kK=U [n] [i+1 ] [j+1 ]; 

lL=U[n][i-l][j+l]; 

mM=U[n] [i-1 ][j-1 ]; 

nN=U[n][i+l][j-l]; 

oO=V[n][i][j+l]; 

pP=V[n][i][j]; 

qQ=V[n][i][j-l]; 

rR=V[n-l][i][j]; 


EE=FEE(Aa,Bb,Cc,Dd,Ee,Ff,Gg,Hh,Ii,Jj,Kk,Ll,Mm,Nn,Oo,Pp,Qq,Rr); 

FF=FFF(Aa,Bb,Cc,Dd,Ee,Ff,Gg,Hh,Ii,Jj,Kk,Ll,Mm,Nn,Oo,Pp,Qq,Rr); 

GG=FGG(aA,bB,cC,dD,eE,fF,gG,hH,iI,jJ,kK,lL,mM,nN,oO,pP,qQ,rR); 

HH=FHH(aA,bB,cC,dD,eE,fF,gG,hH,iI,jJ,kK,lL,mM,nN,oO,pP,qQ,rR); 

/* Calculating the n+1 values */ 
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u [n+1 ] [i] [j ]=(BB *EE-AA*FF)/(BB *CC+A A*DD); 
v[n+l][i][j]=(BB*GG-AA*HH)/(BB*CC+AA*DD); 
U[n+l][i][j]=(DD*EE+CC*FF)/(BB*CC+AA*DD); 
V[n+l][i][j]=(DD*GG+CC*HH)/(BB*CC+AA*DD); 

} 

} 

for (i=0;i<=50;i++) 

{ for (j=0;j<=50;j++) 

{ 

/* Updating the grid for next time level */ 

u[n-l][i][j]=u[ n ][i][j]; 
u[n][i][j]=u[n+l][i]U]; 
v[n-1 ] [i]0]=v[ii] [i] [j]; 
v[n][i][j]= V [n+l][i]|j]; 

U[n-l][i]|J]=U[n][i][j]; 

U [n] [i] [j]=U[n+1 ] [i] [j]; 

V[n-l][i](j]=V[n][i][j]; 

V[n][i][j]=V[n+l][i]|J]; 

} 

} 

for (i=l;i<=49;i++) 

{ 

sol[ind][i]=sqrt(v[2][i][l]*v[2][i][l]+u[2][i][l]*u[2][i][l]); 

flu[ind][i]=sqrt(V[2][i][l]*V[2][i][l]+U[2][i][l]*U[2][i][l]); 

} 


ind = ind + 1; 

) 

for (tt=0;tt<time;tt++) 

{ 

for (i=0;i<=50;i++) 

printf("\n%d\t%d\t%g\t%g",tt,i,sol[tt][i],flu[tt][i]); 


} 

/*for (i=l ;i<=49;i++) 

{ 
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for (j=l;j<=49;j++) 

{ 

printf( "\n%d\t%d\t%g\t%g\t%g\t%g" ,i,j,u[n] [i] [j],v[n] [i] [j], 
U[n][i][j],V[n][i][j]); 


}*/ 

return 0; 

} 


/* Declared Function Definition */ 

float FEE(float Aa,float Bb,float Cc,float Dd.float Ee,float Ff,float Gg, 
float Hh,float Ii,float Jj, float Kk.float LI,float Mm.float Nn, 
float Oo,float Pp,float Qq.float Rr) 

{ 

a = (l-beta+(beta*Kr/Kf)-(Kh/Kr)); 

Q = (beta*( 1 -beta)*Kr-beta*Kb)/a; 

R = beta*beta*Kr/a; 

P = ((l-beta)*((l-beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a + 4*Mu/3; 

b = N*beta*beta/K; 

rhol2 = rhof*beta*(l-alpha); 

rholl = (l-beta)*rhor-rhol2; 

rho22 = rhof*beta-rhol2; 

ans = Q/Mu*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h)) + 

R/Mu*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h)) - 
b*Fw*w*((Dd-Nn)/(2*delT*Mu)) - 

rho22*(P-Mu)/(Mu*rhol2)*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h)) - 
rho22/rho 12*(Aa-4*Bb+Cc+Ee+Ff)/(h*h) - 

rho22*Q/(rhol2*Mu)*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h)) - 
rho22*b*Fw*w*(Dd-Nn)/(rhol2*2*delT*Mu) - 
(Dd-2*Bb)*w*w*(rhol2-((rho22*rholl)/rhol2))/(delT*delT*Mu); 
return ans; 

} 


float FFF(float Aa,float Bb,float Cc,float Dd.float Ee.float Ff,float Gg, 
float Hh,float Ii,float Jj.float Kk,float Ll.float Mm.float Nn, 
float Oo,float Pp.float Qq.float Rr) 

{ 

a = (l-beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q = (beta*( 1 -beta)*Kr-beta*Kb)/a; 

R = beta*beta*Kr/a; 
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P = ((l-beta)*((l-beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a + 4*Mu/3; 

b = N*beta*beta/K; 

rhol2 = rhof*beta*(l -alpha); 

rhol 1 = (l-beta)*rhor-rhol2; 

rho22 = rhof*beta-rhol2; 

ans = Q/Mu*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h)) + 

R/Mu*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h)) - 
b*Fw*w*((Dd-Nn)/(2*delT*Mu)) - 

rho 12*(P-Mu)/(rhol 1 *Mu)*((Aa-2*Bb+Cc)/(h*h)+(Ii-Gg+Hh-Jj)/(4*h*h)) - 
rhol2/rhol 1 *(Aa-4*Bb+Cc+Ee+Ff)/(h*h) - 

rho 12*Q/( rhol l*Mu)*((Kk-2*Ll+Mm)/(h*h)+(Oo-Pp+Qq-Rr)/(4*h*h)) - 
rhol2*b*Fw*w*(Dd-Nn)/(rhol l*2*delT*Mu) - 
(Nn-2*Ll)*w*w*(rho22-((rhol2*rhol2)/rhol l))/(delT*delT*Mu); 
return ans; 


float FGG(float aA,float bB,float cC,float dD,float eE,float flF,float gG, 
float hH,float il,float jJ,float kK,float 1L,float mM,float nN, 
float oO,float pP,float qQ,float rR) 

{ 

a = (1 -beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q = (beta*(l-beta)*Kr-beta*Kb)/a; 

R = beta*beta*Kr/a; 

P = ((l-beta)*((l-beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a + 4* Mu/3; 

b = N*beta*beta/K; 

rho 12 = rhof*beta*(l-alpha); 

rhol 1 = (l-beta)*rhor-rhol2; 

rho22 = rhof*beta-rhol2; 

ans = Q/Mu*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h)) + 
R/Mu*((kK-lL+mM-nN)/(4*h + h)+(oO-2*pP+qQ)/(h*h)) - 
b*Fw*w*(hFl-rR)/(2*delT*Mu) - 

rho22*(P-Mu)/(rhol2*Mu)*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h))- 
rho22/rho 12*(iI-4*fF+jJ+eE+gG)/(h*h) - 

rho22*Q/(rhol2*Mu)*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h)) - 
rho22*b*Fw*w*(hH-rR)/(rhol2*2*delT*Mu) - 
(hH-2*flF)*w*w*(rhol2-((rho22*rhol l)/rhol2))/(delT*delT*Mu); 

} 

float FFlH(float aA,float bB,float cC,float dD,float eE.float flF,float gG, 
float hH,float iflfloatjJ,float kK,float 1L,float mM,float nN, 
float oO,float pP.float qQ,float rR) 

( 
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a = (l-beta+(beta*Kr/Kf)-(Kb/Kr)); 

Q = (beta*(l-beta)*Kr-beta*Kb)/a; 

R = beta*beta*Kr/a; 

P = ((l-beta)*((l-beta)-Kb/Kr)*Kr+(beta*Kr*Kb/Kf))/a + 4*Mu/3; 

b = N*beta*beta/K; 

rhol2 = rhof*beta*(l-alpha); 

rholl = (l-beta)*rhor-rhol2; 

rho22 = rhof*beta-rhol2; 

ans = Q/Mu*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h)) + 

R/Mu*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h)) - 
b*Fw*w*(hH-rR)/(2*delT*Mu) - 

rho 12*(P-Mu)/(rho 11 *Mu)*((aA-bB+cC-dD)/(4*h*h)+(eE-2*fF+gG)/(h*h)) - 
rhol2/rhol l*(iI-4*fF+jJ+eE+gG)/(h*h) - 

rhol2*Q/(rhol l*Mu)*((kK-lL+mM-nN)/(4*h*h)+(oO-2*pP+qQ)/(h*h)) - 
rho 12*b*Fw*w*(hH-rR)/(rho 11 *2*delT*Mu) - 
(rR-2*pP)*w*w*(rho22-((rhol2*rhol2)/rhol l))/(delT*delT*Mu); 
return ans; 

} 
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